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Nonlinear Particle Acceleration in Relativistic Shocks 
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ABSTRACT 

Monte Carlo techniques are used to model nonlinear particle acceleration in paral- 
Ph. lei collisionless shocks of various speeds, including mildly relativistic ones. When the 
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acceleration is efficient, the backreaction of accelerated particles modifies the shock 



structure and causes the compression ratio, r, to increase above test-particle values. 
Modified shocks with Lorentz factors, 70 < 3, can have compression ratios considerably 
greater than 3 and the momentum distribution of energetic particles no longer follows 
\q • a power law relation. These results may be important for the interpretation of gamma- 

ray bursts if mildly relativistic internal and/or afterglow shocks play an important role 
. accelerating particles that produce the observed radiation. For 70 > 10, r approaches 

3 and the so-called 'universal' test-particle result of N(E) oc E~ 2 - 3 is obtained for 
sufficiently energetic particles. In all cases, the absolute normalization of the particle 
distribution follows directly from our model assumptions and is explicitly determined. 

Subject headings: Cosmic rays — acceleration of particles — relativistic shock waves — 
& . gamma-ray bursts; PACS: 52.60, 96.40 



1. Introduction 

Most collisionless shocks in astrophysics are nonrelativistic, i.e., the flow speed of the unshocked 
plasma in the reference frame at rest with the shock, no, is much less than the speed of light, c. 
Particle acceleration in such shocks has been studied extensively in both the linear and nonlinear 
regimes, and we refer the reader to reviews which describe the basic features of the shocks, the 
energetic particles they produce, and the many applications (e.g., Drury 1983; Blandford & Eichler 
1987; Jones &: Ellison 1991). Relativistic shocks, where the flow speed Lorentz factor 70 = [1 — 
(«o/c) 2 ] -1 / 2 is greater than a few, are likely to be much less common than nonrelativistic ones, but 
may occur in extreme objects such as pulsar winds, hot spots in radio galaxies, and gamma-ray 
bursts (GRBs). Largely motivated by the application to GRBs, relativistic shocks have recently 
received considerable attention by a number of researchers (e.g., Bednarz & Ostrowski 1996; Kirk 
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et al. 2000; Achterberg et al. 2001). However, except for some preliminary work done over a 
decade ago (Schneider k Kirk 1987; Ellison 1991a, b), current descriptions of relativistic shocks 
undergoing first-order Fermi acceleration are test particle approximations that do not include the 
backreaction of the accelerated particles on the shock structure. This may be a serious limitation of 
relativistic shock theory in applications, such as GRBs, where high particle acceleration efficiencies 
are often assumed. Here, we present results for nonlinear acceleration where the backreaction of 
the accelerated population on the relativistic shock structure is included self-consistently. 

In collisionless shocks, charged particles gain energy by scattering back and forth between the 
converging upstream and downstream plasmas. This basic physical process, called diffusive or first- 
order Fermi shock acceleration, is the same in relativistic and nonrelativistic shocks. Differences 
in the mathematical description and outcome of the process occur, however, because energetic 
particle distributions are nearly isotropic in the shock reference frame in nonrelativistic shocks 
(where v S> uq; v is the individual particle speed), but are highly anisotropic in relativistic shocks 
(since v ~ no ~ c) (e.g., Peacock 1981). The description of particle diffusion and energy gain is far 
more difficult when 70 3> 1 because the diffusion approximation, which requires nearly isotropic 
distribution functions, cannot be made. Because of this, Monte Carlo simulations, where particle 
scattering and transport are treated explicitly, and which, in effect, solve the Boltzmann equation 
with collective scattering (e.g., Ellison k Eichler 1984; Kirk k Schneider 1987b; Ellison, Jones, 
k Reynolds 1990; Ellison k Reynolds 1991; Ostrowski 1991; Bednarz k Ostrowski 1996), offer 
advantages over analytic methods. This is true in the test-particle approximation, where analytic 
results exist, but is even more important for nonlinear relativistic shocks. 

In nonrelativistic shocks, for v 2> uq, a diffusion-convection equation can be solved directly for 
infinite, plane shocks (e.g., Axford, Leer, k Skadron 1977; Blandford k Ostriker 1978), yielding the 
well-known result 

f{p)dPpctp~ a d' i p with a = 3r/(r — 1) , (1) 

where r is the shock compression ratio, p is the momentum, and f(p) d 3 p is the number density of 
particles in d?p. Eq. (1) is a steady-state, test-particle result with an undetermined normalization, 
but the spectral index, a, in this limit is independent of the shock speed, uo, or any details of the 
scattering process as long as there is enough scattering to maintain isotropy in the local frame. To 
obtain an absolute injection efficiency, or to self-consistently describe the nonlinear backreaction 
of accelerated particles on the shock structure (at least when the seed particles for acceleration 
are not fully relativistic to begin with), techniques which do not require « > «o must be used. 
Furthermore, for particles that do not obey v ^> uq additional assumptions must be made for how 
these particles interact with the background magnetic waves and/or turbulence, i.e., the so-called 
"injection problem" must be considered (see, for example, Jones k Ellison 1991; Malkov 1998). The 
Monte Carlo techniques we describe here make the simple assumption that all particles, regardless 
of energy, interact in the same way, i.e., all particles scatter elastically and isotropically in the 
local plasma frame with a mean free path proportional to their gyroradius. These techniques and 
assumptions have been used to calculate nonlinear effects in nonrelativistic collisionless shocks for a 
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number of years with good success comparing model results to spacecraft observations (e.g., Ellison 
& Eichler 1984; Ellison, Mobius, & Paschmann 1990; Ellison, Jones, & Baring 1999). 

Early work on relativistic shocks was mostly analytical in the test particle approximation (e.g., 
Blandford & McKee 1976; Peacock 1981; Kirk & Schneider 1987a; Heavens & Drury 1988), although 
the analytical work of Schneider & Kirk (1987) explored modified shocks. Test-particle Monte Carlo 
techniques for relativistic shocks were developed by Kirk & Schneider (1987b) and Ellison, Jones, 
& Reynolds (1990) for parallel, steady-state shocks, i.e., those where the shock normal is parallel to 
the upstream magnetic field, and extended to include oblique magnetic fields by Ostrowski (1991). 
Some preliminary work on modified relativistic shocks using Monte Carlo techniques was done by 
Ellison (1991a,b). 

The most important results from the theory of test-particle acceleration in ultrarelativistic 
shocks are: (i) regardless of the state of the unshocked plasma, some particles can pick up large 
amounts of energy AE ~ 7q mc 2 in their first shock crossing cycle (Vietri 1995), but will receive 
much smaller energy boosts (<Ej / 'Ei><^ 2) for subsequent crossing cycles (e.g., Gallant h Achter- 
berg 1999; Achterberg et al. 2001) 2 ; (ii) the shock compression ratio, defined as r = uq/u2, tends 
to 3 as uo — > c (e.g., Peacock 1981; Kirk 1988), where U2 is the flow speed of the shocked plasma 
measured in the shock frame; 3 and (iii) a so-called 'universal' spectral index, a ~ 4.2 — 4.3 (in 
Eq. 1) exists in the limits of 70 S> 1 and 56 <C 1, where 59 is the change in direction a particle's mo- 
mentum vector makes at each pitch angle scattering (e.g., Bednarz & Ostrowski 1998; Achterberg 
et al. 2001). 

We find that these results are modified in mildly relativistic shocks, even in the test-particle 
approximation, and in fully relativistic shocks (at least for 70 < 10) when the backreaction of 
the accelerated particles is treated self-consistently, which causes the shock to smooth and the 
compression ratio to change from test-particle values. In mildly relativistic shocks, f(p) remains a 
power law in the test-particle approximation but both r and a depend on the shock Lorentz factor, 
7o- When efficient particle acceleration occurs in mildly relativistic shocks (i.e., 70 < 3), large 
increases in r can result and a power law is no longer a good approximation to the spectral shape. 
In these cases, we determine the compression ratio by balancing the momentum and energy fluxes 
across the shock with the Monte Carlo simulation. For larger Lorentz factors, accelerated particles 
smooth the shock structure just as they do in slower shocks, but r approaches 3 as 70 increases. In 
general, efficient particle acceleration results in spectra very different from the so-called 'universal' 
power law found in the test-particle approximation unless 70 > 10. 



2 Ei (Ef) is the particle energy at the start (end) of an upstream to downstream to upstream (or a downstream 
to upstream to downstream) shock crossing cycle. 

3 Note that the density ratio across the relativistic shock P2/P0 =fc uo/112, in contrast with nonrelativistic shocks, 
because the Lorentz factors associated with the relativistic flows modify the particle flux jump condition. Here and 
elsewhere we use the subscript (2) to indicate far upstream (downstream) values. 
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2. Monte Carlo Model 

The techniques we use are essentially identical to those described in Ellison, Baring, & Jones 
(1996) and Ellison, Jones, & Baring (1999). The differences are that the code has been made fully 
relativistic and only results for parallel shocks with pitch-angle diffusion are presented here. The 
code is steady-state, includes a uniform magnetic field, and moves particles in helical orbits. We 
assume the Alfven Mach number is large, i.e., we neglect any effects from Alfven wave heating 
in the upstream precursor. This also means we neglect the second-order acceleration of particles 
scattering between oppositely propagating Alfven waves. Such an effect in relativistic plasmas with 
strong magnetic fields is proposed for nonlinear particle acceleration in GRBs by Pelletier (1999) 
(see also Pelletier k Marcowith 1998). 

The pitch angle diffusion is performed as described in (Ellison, Jones, & Reynolds 1990). That 
is, after a small increment of time, St, a particles' momentum vector, p, undergoes a small change in 
direction, 59. If the particle originally had a pitch angle, 9 (measured relative to the shock normal 
direction), it will have a new pitch angle 9' such that 

cos 9' = cos 6 cos 59 + \/l — cos 2 9 sin 59 cos 4> , (2) 

where 4> is the azimuth angle measured with respect to the original momentum direction. All angles 
are measured in the local plasma frame. If 59 is chosen randomly from a uniform distribution 
between and 59 max and <\> is chosen from a uniform distribution between and 2tt, the tip of the 
momentum vector will perform a random walk on a sphere of radius p. As shown by Ellison, Jones, 
& Reynolds (1990), the angle 59 ma , x is determined by 

59 max = (6 5t/t c ) 1/2 = {Utt/N) 1 ' 2 , (3) 

where N = T g /5t 3> 1 is the number of gyro-segments, 5t, dividing a gyro-period T g = 2nr g /v. The 
time t c is a "turn around" time defined as t c = X/v, where A is the particle mean free path. The 
mean free path is taken to be proportional to the gyroradius r g = pc/{QeB) (e is the electronic 
charge, Q is the ionic charge number, and B is the local uniform magnetic field), i.e., A = rjr g , 
where r\ determines the strength of scattering. In all of the examples given here we set rj = 1, or, 
in other words, the strong scattering Bohm limit is assumed. 

For a downstream particle to return upstream, it's velocity vector must be directed within 
a cone with opening angle 92 such that jt^cos^l > U2, where V2 and 92 are measured in the 
downstream frame and 92 = 0° is in the — x-direction, i.e., along the shock normal direction. For 
fully relativistic shocks with v 2 — c and 112 — c/3, cos #2 <^ 1/3 for a downstream particle to cross 
the shock into the upstream region. When the particle enters the upstream region it must satisfy 
essentially the same constraint, i.e., |t>ocos#o| > uq, where now vq and 9q are measured in the 
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upstream frame. 4 Since both the particle and shock have high Lorentz factors, we can write 

where j v = [1 — (t>o/c) 2 ] -1 / 2 is the particle Lorentz factor. Since cos#o — 1 — 9q/2 for small 9q, we 
have 

(5) 

For ultrarelativistic particles with j v » 70, #0 — I/70 ( e -g-) Gallant & Achterberg 1999), but 9q 
can be much smaller for mildly relativistic particles. 

In order to re-cross into the downstream region, particles must scatter out of the upstream cone 
defined by #0 an d Achterberg et al. (2001) show that most particles are only able to change the 
angle they make with the upstream directed shock normal by \59\ ~ #0 before being sweep back 
downstream, making the distribution of shock crossing particles highly anisotropic. Therefore, if 
the shock Lorentz factor 70 3> 1, a larger fraction of particles re-cross the shock into the downstream 
region with highly oblique angles (as measured in the shock frame) compared to lower speed shocks 
(see Fig. 2 discussed below). Particles crossing at such oblique angles receive smaller energy gains 
than would be the case for an isotropic pitch angle distribution and Achterberg et al. (2001) go on 
to show that <Ef/Ei>~ 2 for a shock crossing cycle (after the first one). 

Considering these constraints, we require 59 max < 9q, or 

N > N min = 12vr 7o 2 , (6) 

and find (as shown in Fig. 1) that the power law spectral index, a, asymptotically approaches a 
maximum value as N is increased. If N is less than the value required for convergence (and the gyro- 
segments are too large), the distribution will be flatter than produced with the convergent value 
of N because more particles are able to cross from upstream to downstream with 9^ <C 90° and 
receive unrealistically large energy boosts. This effect has long been known from the comparison of 
pitch-angle diffusion to large-angle scattering in relativistic shocks (e.g., Kirk & Schneider 1987b; 
Ellison, Jones, & Reynolds 1990). For all of the examples reported on here, ./V is chosen large 
enough so it makes no difference if 59 is chosen uniformly between 0° and <50 max or if cos 59 is 
chosen uniformly between cos <56> max and 1. Fig. 1 shows how our results depend on iV for shock 
speeds ranging from fully nonrelativistic to fully relativistic. In all cases, as iV is increased the 
spectral index approaches a maximum and for 70 > 7 we obtain the well known result a = 4.2-4.3. 
The fact that the computation time for the Monte Carlo simulation scales as ./V and N oc 7 2 places 
limits on modeling ultrarelativistic shocks with this technique. 



4 In the test-particle approximation, uo is just the shock speed. In nonlinear shocks, the flow speed just upstream 
from the subshock at x = will be less than the far upstream shock speed, Mo, as measured in the shock reference 
frame. 
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Fig. 1. — Power law spectral index, a, versus number of gyro-segments, N, for various test-particle 
shocks as labeled. In all cases, as N is increased the spectral index converges and we obtain the 
known results of a ~ 4 for nonrelativistic shocks with r = 4, and a ~ 4.23 for fully relativistic 
shocks with r = 3. The parameter iV m i n is defined in Eq. (6). 
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Fig. 2. — Distribution of the cosine of the pitch angle, i.e., p x /pt, for particles crossing x = 0. The 
solid and dashed curves in the top panel are for a shock with Lorentz factor 70 = 2. The dotted 
curve in the top panel is for a nonrelativistic shock with speed no = 5000 km s _1 . The bottom 
panel shows curves for a fully relativistic shock with 70 = 10. In all cases, the curves are normalized 
so that the areas under them is 1 and the x-component of momentum, p x , is positive when directed 
downstream. The nonlinear (NL) and unmodified (UM) shock results are labeled. 
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In Fig. 2 we compare pitch-angle distributions (measured in the shock reference frame) of 
particles crossing the shock. The curves are normalized such that the area under each curve equals 
one and the x-component of particle momentum in the shock frame, p x , is positive when directed 
downstream to the right (see Fig. 8 for the shock geometry). Here, pt is the magnitude of the 
total particle momentum also measured in the shock frame. In the top panel, we compare an 
unmodified (UM) and nonlinear (NL) mildly relativistic shock (70 = 2) with a nonrelativistic one 
(uo = 5000 km s -1 ). Particles crossing the 70 = 2 shock are highly anisotropic with p x /pt strongly 
peaked near ~ 0.35. In the nonrelativistic shock, the particles are nearly isotropic except for a 
slight flux- weighting effect. There is little difference in the distributions between the UM and NL 
shocks. In the bottom panel, we show the pitch-angle distributions for UM and NL shocks with 
70 = 10. While the distributions are somewhat more sharply peaked, they are quite similar to 
those for 70 = 2 and show relatively small variations between the UM and NL shocks. 

The main difference between our present code and an earlier code used by Ellison, Jones, & 
Reynolds (1990) to model test-particle relativistic shocks is that the previous code used a guiding 
center approximation with an emphasis on large-angle scattering rather than the more explicit 
orbit calculation of pitch-angle diffusion used here. Other than the far greater range in 70 and the 
nonlinear results we now present, the work of Ellison, Jones, & Reynolds (1990) is consistent with 
the work presented here. 

The particle transport is performed as follows. Particles of some momentum, p p f, (measured 
in the local plasma frame) are injected far upstream from the shock and pitch-angle diffuse and 
convect until they cross a grid zone boundary, i.e., a dividing plane in the simulation between 
regions with different bulk flow speeds (see Ellison, Baring, & Jones 1996, for a full discussion). For 
unmodified shocks there is only one grid zone boundary which divides the upstream and downstream 
regions, but for nonlinear shocks the bulk flow speed changes in small steps, each separated by a 
boundary, from uq far upstream to ui downstream. In the unmodified case, the shock thickness 
is essentially zero (i.e., shorter than the distance a particle diffuses in St) but in the nonlinear, 
modified case, the shock precursor extends over the entire region of varying bulk flow speeds and 
a small scale "subshock" (at position x = in our simulation) exists where most of the entropy 
production occurs. When a particle crosses a grid zone boundary, p p f is transformed to the new 
local frame moving with a new speed relative to the subshock and the particle continues to scatter 
and convect. Each particle is followed until it leaves the system in one of three ways. It can convect 
far downstream and not return to the subshock, it can obtain a momentum greater than some p max 
and be removed, or it can diffuse far enough upstream to cross an upstream free escape boundary 
(FEB) and be removed. Both p max and the position of the FEB are free parameters in our model 
(see Berezhko & Volk 1997, for a discussion of the self-consistent determination of the maximum 
particle energy in nonrelativistic supernova remnant shocks). For our nonlinear calculations, we 
iterate the shock structure and compression ratio until the number, momentum, and energy fluxes 
are conserved across the shock. This procedure has been detailed many times for nonrelativistic 
shocks (see Ellison, Baring, & Jones 1996, and references therein) and the modifications required 
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for relativistic shocks were given in Ellison & Reynolds (1991). 

To avoid excessive computation, we use a probability of return calculation as described in detail 
in Ellison, Jones, & Reynolds (1990). That is, we use the standard expression obtained by Peacock 
(1981), 

P R =(^)\ (7) 

to determine the probability, Pr, that a particle, having crossed a particular point in the uniform 
downstream flow, will return back across that point. Eq. (7) is fully relativistic and independent of 
the diffusive properties of the particles as long as they are isotropic in the u 2 frame. We ensure this 
isotropy by only applying Eq. (7) once a particle has diffused several mean free paths downstream 
from the shock. For ultrarelativistic shocks with v ~ uq ~ c and r ~ 3, Eq. (7) gives Pr ~ 0.25. 

For clarity, we note that this is not the same probability, V Te t, that is used by Achterberg et al. 
(2001) to determine the test-particle power law index, i.e., 

ln(l/P ret ) 

+ In <E f /Ei> ' [ ' 

where N(E) oc E~ s and, for fully relativistic particles, s = a — 2. The quantity <Ef/Ei> is 
the average energy ratio for a particle undergoing a shock crossing cycle. Although not explicitly 
stated in Achterberg et al. (2001), it is clear from the context that V Te t is calculated just behind 
the shock where the particle distribution is highly anisotropic. In this case, particles that have 
just crossed from upstream to downstream will be more likely to recross back into the upstream 
region than indicated by Eq. (7) because their pitch angles are more likely to be highly oblique 
relative to the shock normal than in the isotropic distributions further downstream (compare the 
solid or dashed curves to the dotted curve in the top panel of Fig. 2 discussed below). For 70 = 10, 
Achterberg et al. (2001) find V rct = 0.435 ± 0.005 and <E f /Ei>= 1.97 ± 0.01 giving the standard 
result s = 2.230±0.012. As shown in Fig. 1, our unmodified results are consistent with this spectral 
index for 70 > 7. 



3. Test-Particle Results 

In Fig. 3 we show particle distributions for unmodified shocks with speeds ranging from fully 
nonrelativistic (no = 5000 km s -1 ) to mildly relativistic {(5q = uq/c = 0.5) to fully relativistic 
(70 = 10). The nonrelativistic distribution matches the standard test-particle Fermi result of a = 4 
for r = 4 and the fully relativistic result is consistent with the well-known limit of a — ► 4.2 — 4.3 
as 70 — > 00 for r = 3. In the trans-relativistic regime, both the compression ratio and the spectral 
index vary with uq. For the uq = 0.5c distribution shown in Fig. 3, we have determined the 
compression ratio by balancing the mass, momentum, and energy fluxes across the shock in the 
test-particle limit, i.e., by ignoring any effects from the accelerated particles. This technique is 
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Fig. 3. — Particle spectra, p 4 f(p), versus momentum for various unmodified shocks with speeds 
as indicated. The test-particle compression ratios, r, and spectral indices, a, are noted. The far 
upstream plasmas in the 5000 km s _1 , and (3q = uq/c = 0.5 shocks are thermal at 10 6 K. The 
7o = 10 shock has a far upstream plasma which is a ^-function at 1 MeV. All spectra are calculated 
at the shock, in the shock frame, and the relative normalization is arbitrary. 



described in detail in Ellison & Reynolds (1991). We find r = 3.8 ± 0.1 and the spectral index is 
a ~ 4.01, slightly flatter than 3r/(r — 1) ~ 4.07, the nonrelativistic result. 
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Fig. 4. — The solid line is the compression ratio, r, and the dashed line is the spectral index, a, for 
unmodified (i.e., test-particle) shocks. The solid dots show r for shocks undergoing efficient particle 
acceleration. In all cases, r is determined for each /3o7o by balancing the momentum and energy 
fluxes across the shock. The maximum cutoff momentum, p maX ) is unimportant for the test-particle 
shocks, but is set to 10 4 /?o7o "ipC for the nonlinear shocks. 

Fig. 4 shows the compression ratio as a function of /3o7o (solid curve), still ignoring the effects 
of accelerated particles. The compression ratio is determined self-consistently by balancing the 
momentum and energy fluxes across the shock with no restriction on the shocked or unshocked 
adiabatic index. As expected, r decreases smoothly from 4 for fully nonrelativistic shocks to ~ 3 
for fully relativistic shocks. The power law index, a, is also shown (dashed curve) and this varies 
slowly from a = 4 to a ~ 4.23 between the two extremes. For comparison, we show r (solid dots) 
for shocks undergoing efficient particle acceleration. For these points, all shock parameters except 
/3o7o an d Pm&x are kept constant. For the nonlinear shocks, the maximum cutoff momentum is set 
to p ma x = 10 /?o7o ?™pC in all cases. 5 We discuss the nonlinear results in detail below, but here only 



'The value of Pmax can have a large influence on the shock characteristics at low /3o7o because r is large enough 
(and a small enough) that particles escaping at p max are dynamically important. When /3o7o becomes large enough 
so that r drops below ~ 4, p max is no longer an important parameter for the shock structure. 
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cmphasize that r in nonlinear shocks will be larger than the test-particle value for 70 < 10. The 
test-particle results shown in Fig. 4 are in close agreement with recent analytic results of Kirk et 
al. (2000) and Gallant (2002). To our knowledge, no analytic results exist for nonlinear, relativistic 
shocks. 
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Fig. 5. — Average ratios of final (f) to initial (i) momentum for downstream to upstream to 
downstream and upstream to downstream to upstream shock crossing cycles. The histograms are for 
a test-particle shock with 70 = 10 and r = 3. Note the large momentum gain (<p f lp{> u -d-u= 115) 
in the first shock crossing cycle. Note also that this particular value depends on our choice of the 
injection momentum. 

In Fig. 5 we show the average ratios of momenta (measured in the local plasma frame) for 
particles executing upstream to downstream to upstream cycles across the shock, <p//pi> u -d-u, 
and downstream to upstream to downstream cycles, <p f /pi> d-u-d- These results are for 70 = 10 
(r = 3) and show a slight momentum dependence at low momenta but converge to <Pf /pi>= 
2.2 ± 0.1 at high momenta. This value is close to <Ef /Ei>= 1.97 ± 0.01 reported by Achterberg 
et al. (2001) for 70 = 10. As mentioned above, in the first upstream to downstream to upstream 
cycle, particles achieve a large boost in momentum as indicated in the figure. 

The difference between our value of <Pf/pi> at large pi and that of Achterberg et al. (2001) 
is greater than the uncertainties and probably stems from the different assumptions made in the 
simulations. As discussed above, <Pf/pi> depends critically on the average angle a particle makes 



1 1 II 1 
- 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

1 

^<P f /Pi> u _ d _ u =115 j 




| Dw— Up— Dw 

H 'L_ 
i i 




! I 

L 1 

\Up-Dw-Up 


1 1 1 1 1 


7 =10, r=3 " 



-13- 



when crossing the shock. While the large majority of particles in our simulations gain energy 
by crossing from downstream to upstream and then immediately (within a few <5i's) re-crossing 
back into the downstream region at oblique angles, a few manage to diffuse farther upstream (see 
Fig. 7 below). When these particles re-cross the shock into the downstream region, they can do so 
at flatter angles and receive larger energy gains. We speculate that differences in how these few 
particles are treated in the simulations might produce the differences in <Pf/Pi>- 

As an illustration of how particles interact with nonrelativistic and relativistic unmodified 
shocks, we shown trajectories for two individual particles in Figs. 6 and 7. The lower panel in each 
figure shows a trace of the particle trajectory and the upper panels show the particle momentum, 
always measured in the local plasma frame. For the nonrelativistic shock (Fig. 6), the speed of the 
particle is far greater than the shock speed and it diffuses easily on both sides of the shock. When 
it crosses x = 0, it does so nearly isotropically (except flux weighting makes crossings with flat 
trajectories slightly more likely) and essentially always gains momentum. The momentum gain in 
a single shock crossing is small, but a particle can stay in the system for many crossings. 

When the shock speed, uq, is close to c, the particle will be convected downstream much more 
rapidly than in nonrelativistic shocks and few particles will be able to cross the shock many times. 
However, downstream particles that do manage to cross the shock into the upstream region do so 
with much flatter trajectories, as discussed above, and can receive large momentum boosts in a 
single shock crossing due to the shock's Lorentz factor (note the logarithmic scale in the top panel 
of Fig. 7). In a typical shock crossing cycle, downstream particles gain momentum when they cross 
into the upstream region (see positions labeled a and b in Fig. 7), lose momentum when they cross 
back downstream because they cross with oblique pitch angles, but end up with a net momentum 
boost. However, as shown by the position labeled c, it is possible for a particle to diffuse farther 
upstream before being convected back to the shock. In this case, it can cross the shock with a 
flat trajectory and gain momentum upon entering the downstream region. If the acceleration is 
efficient, the few particles that diffuse far upstream carry enough pressure to produce the shock 
smoothing we discuss next. 



4. Non-Linear Results 

As noted by Achterberg et al. (2001), the large energy boost particles receive in their initial 
crossing of the shock provides a natural injection process for further acceleration and suggests 
that relativistic shocks may be efficient accelerators. However, just as with nonrelativistic shocks, 
efficient acceleration limits the use of test-particle approximations and requires that the nonlinear 
backreaction of the accelerated particles be treated self-consistently (e.g., Jones & Ellison 1991). 
These nonlinear effects will result in a smoothing of the shock and a change in the overall shock 
compression ratio, just as they do in nonrelativistic shocks. The differences between test-particle 
and nonlinear results, for the parameter ranges we have investigated, are large enough to produce 
spectra noticeably different from the often quoted N(E) oc E~ 2 - 3 and to influence applications to 
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Fig. 6. — Particle trajectory (lower panel) and momentum (upper panel) in an unmodified shock of 
speed uq = 5000 km s _1 . The momentum is calculated in the local plasma frame, either upstream 
or downstream from the shock. 
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Fig. 7. — Trajectory and momentum for a particle in an unmodified shock with speed uq = 0.9c. 
Note that the horizontal axis is split at 1200 scatterings. 
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GRB models where high particle acceleration efficiencies are assumed. We illustrate this with two 
examples, one mildly relativistic (70 = 1.4) and one more fully relativistic (70 = 10). The far 
upstream conditions have relatively little influence on our results as long as 7 P <C 70, where j p is 
the plasma frame Lorentz factor of the far upstream, injected particles. For concreteness, in both 
examples we take the far upstream plasma to be a thermal distribution of protons at a temperature 
of 10 6 K. 



4.1. Mildly relativistic shock: 70 = 1.4 
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Fig. 8. — Unmodified (UM) and nonlinear (NL) shock profiles, i.e., j u (x)u(x), and momentum and 
energy fluxes versus position, x. All quantities are scaled to far upstream values and in all panels 
the solid curves are results from unmodified shocks and the dashed curves are nonlinear results. 
The NL momentum and energy fluxes are 3 to 4% below the far upstream values because particles 
escape at a maximum momentum, p max = 10 4 /?o7o r "-pC ~ 10 13 eV/c. We have used N/N m \ D ~ 10 
in both cases. 
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Fig. 8 shows unmodified and nonlinear shock structures for 70 = 1.4. The top panel shows 
7«(x) u(x), where j u ( x ) = {1 ~~ [ti(x)/c] 2 } -1 / 2 , the middle panel is momentum flux, and the bottom 
panel is the energy flux, all scaled to far upstream values. All curves are plotted versus x, where 
x = is the position of the sharp subshock. A logarithmic scale is used for x < — 10r 9i o and a 
linear scale is used for x > — 10r S) o, where r 9j0 = m p uo/(eB). 6 In each panel, the solid curve is 
from an unmodified shock with r ~ 3.6 (see Fig. 4), while the dashed curve is the momentum and 
energy flux conserving result. 

For our pitch angle diffusion model, where particles interact elastically and isotropically in the 
local frame according to Eq. (2) and the discussion following it, particles are accelerated efficiently 
enough at the unmodified shock that the momentum and energy fluxes are not conserved and rise 
well above the allowed far upstream values. In order to conserve these fluxes, the shock structure 
must be smoothed and the overall compression ratio increased above the test-particle value. Our 
computational scheme calculates this compression ratio and flux conserving profile and the result 
is the dashed curve in the top panel with the corresponding momentum and energy fluxes in the 
middle and bottom panels. 

The source of the non-conservation of momentum and energy is the efficient acceleration of 
particles by the sharp flow speed discontinuity. While the actual injection and acceleration effi- 
ciency depends on our particular pitch angle diffusion model, once our scattering assumptions are 
made, the kinematics determine the injection and acceleration of the particles without additional 
parameters. 7 Of course it would have been possible to make assumptions which resulted in an ac- 
celeration efficiency low enough that momentum and energy are approximately conserved without 
significantly smoothing the shock structure or changing the compression ratio from the test-particle 
value. For example, we could have only allowed shocked particles above some Lorentz factor 7; n j 
to recross into the upstream region or arbitrarily restricted the number of particles that recrossed 
into the upstream to a small fraction, /j n j, of all downstream particles. By making /; n j low enough 
or 7i n j high enough we could make the efficiency as low as we wanted. However, there are at least 
three reasons for not making such an assumption. The first is that restricting the acceleration effi- 
ciency requires additional parameters (i.e., 7i n j and/or /i n j) to those needed to describe pitch angle 
diffusion. 8 The second is that models with inefficient acceleration will not help explain GRBs (or 
other objects) that require high efficiencies. If relativistic shocks are inefficient accelerators they 
are not very interesting. If they are efficient, they will have a qualitative resemblance to the results 



6 In these parallel shock simulations where we ignore Alfven wave production, the magnetic field, B, has no other 
effect than setting arbitrary length and time scales. 

7 Other input parameters, such as the Mach number and p m ax or the position of the FEB, influence the accel- 
eration efficiency, but these are "environmental" parameters rather than parameters needed to describe the plasma 
interactions. 

Alternatively, a far more complex model of the plasma interactions can be postulated than done here, inevitably 
requiring additional parameters (e.g., Malkov 1998). 
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we show even if, as is likely, the actual plasma processes are far more complex than the simple 
model we use. A third, perhaps less compelling, reason is that identical scattering assumptions 
as used here have been used for some time in nonrelativistic shocks and shown to match both 
spacecraft observations (e.g., Ellison, Mobius, &; Paschmann 1990) and hybrid plasma simulations 
(e.g., Ellison et al. 1993) of collisionless shocks. 

The increase in compression ratio from r ~ 3.6 to ~ 4.5 shown in Fig. 8 comes about, in 
part, because the particles escaping at p max carry away momentum and energy fluxes which make 
the shocked plasma more compressible. This effect is countered to some degree by the fact that 
the self-consistent shock produces a downstream distribution with a smaller fraction of relativistic 
particles than the unmodified shock so that the downstream adiabatic index T > 4/3. This tends 
to produce a smaller compression ratio. The escaping fluxes show up as a lowering of the dashed 
curves below the far upstream values, as shown in the bottom two panels of Fig. 8, and amount 
to about 3% of the far upstream values for both momentum and energy. Including the escaping 
fluxes, the momentum and energy fluxes are conserved to within a few percent of the far upstream 
values. While the changes seen here are similar to those seen and discussed for many years in 
efficient, nonrelativistic shock acceleration (see Berezhko & Ellison 1999, and references therein), it 
must be noted that there are no known analytic expressions relating escaping fluxes and T to r in 
this trans-relativistic regime. One obvious difference is that for nonrelativistic shocks, the escaping 
momentum flux is generally much less than the escaping energy flux (when both are measured as 
fractions of incoming flux) since p e v1/{povfy S> p e v1/{pou^) when v e 2> no (see Ellison 1985). Here, 
v e is the velocity of the escaping particle. In relativistic shocks, v e ~ uq ~ c so the escaping fluxes 
are about equal as shown in Fig. 8. 

In Fig. 9 we plot p A f{p) for our 70 = 1.4 shock. The solid curve is from an unmodified (UM) 
shock and the dashed curve is the nonlinear (NL) result. The shock smoothing and increase in r 
produce substantial differences in the spectra even though both shocks have exactly the same input 
conditions, (i) The overall normalization of the NL spectrum is less, reflecting the conservation of 
energy flux, (ii) The NL result has the distinctive concave curvature seen in nonrelativistic shocks 
stemming from the fact that higher momentum particles have a longer upstream diffusion length 
and get accelerated more efficiently than lower momentum particles in the smooth shock, (iii) The 
slope at the highest momentum in the NL spectrum reflects the overall compression ratio and is 
flatter than the TP spectrum because r is greater, (iv) The "thermal" peak is shifted to lower 
momentum in the NL result and contains a larger fraction of mildly relativistic particles than in 
the UM result, i.e., T ~ 1.41 for the NL shock compared to T ~ 1.36 for the UM shock. 



4.2. Fully relativistic, nonlinear shock: 70 = 10 

Fig. 10 shows results for 70 = 10. In all panels, the solid curves are for the unmodified shock, 
while the dashed curves are for the flux conserving smoothed shock, both with identical input 
parameters. As before, the acceleration is truncated with a p max = W A /3o^om p c ~ 9.3 x 10 13 eV/c. 
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Fig. 9. — Particle distributions, p 4 f(p), for the shocks shown in Fig. 8. The nonlinear spectrum 
(dashed curve) shows the distinctive concave shape seen in efficient nonrelativistic shock accelera- 
tion, and has a greater fraction of low momentum particles than the spectrum from the unmodified 
shock. As in Fig. 3, the spectra are calculated at the shock in the shock frame and truncated with 
a Pmax- Unlike Fig. 3, the normalization here shows the actual acceleration efficiency. 
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In Fig. 11, u(x) and 7„(x) are plotted separately. Even though the length-scale of the shock 
smoothing is only a few r 9j o and the change in r = 3.03 ± 0.01 is small, they bring the momentum 
and energy fluxes into balance from values a factor of five too large in the unmodified shock. The 
~ 1% difference between r = 3.03 ±0.01 and the canonical value of 3 is significant; a self-consistent 
solution does not exist outside of this range. 
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Fig. 11. — The top panel is the flow speed at x normalized to the far upstream shock speed, no 
versus x, for the 70 = 10 shocks shown in Fig. 10. The bottom panel is the flow Lorentz factor, 
7u(ar), versus x. In both panels, the solid curves are the unmodified shock results with r = 3 and 
the dashed curves are the nonlinear results with r ~ 3.03 ± 0.01. 

As discussed above, both the shock structure and the overall compression ratio may be modified 
when particle acceleration is efficient. For mildly relativistic shocks, such as our 70 = 1.4 example, 
a unique solution can be determined directly from the conserved momentum and energy fluxes. For 
larger 70 's however, the fluxes are less sensitive to changes in the structure and compression ratio 
and an additional constraint is needed to obtain a unique solution. This constraint is provided by 
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the relativistic jump conditions for momentum and energy, i.e., 

2 2 

ll w ^ + p o = J + P2 ; (9) 

lloWQUO = 7«2^2«2 , (10) 

where w = e + P is the enthalpy density, e is the total energy density, P is the pressure, and 
we assume escaping fluxes are negligible. The energy density and pressure are related through a 
combination of the adiabatic equation of state and the conservation of energy, i.e., 

P = (T-l)(e-pc 2 ) , (11) 

where pc 2 is the rest mass energy density and T is, in special cases, the ratio of specific heats (e.g., 
Ellison & Reynolds 1991). Dividing Eq. (9) by Eq. (10) yields 

^ + ^ = ^ + ^ (12) 

c 2 ll {e + P )u c 2 7« 2 (e2 + p 2)«2 ' 



or, in terms of beta's, 



A + ,_5. ^iz£_ A + f_*Jiz£. (13) 



The second term on the left hand side is small compared to (5q for unshocked upstream particle 
temperatures less than 10 9 K and upstream particle densities less than 100 cm -3 . Neglecting this 
term, Eq. (13) becomes 

ft _ A + (_?Uiz/«, (14) 



,e 2 + P 2 / fa 
or, 

e 2 /? 2 2 -/3 (e 2 + P 2 )/3 2 + P 2 = , (15) 



which has the shock solution 



/3o(e 2 + P 2 ) - V^o^ + P 2 ) - 4e 2 P 2 



(16) 



F° r 7o ^ 10, /?o — 1 an d r — e 2 /P 2 . Furthermore, if ultrarelativistic downstream particle speeds 
can be assumed so that e 2 3> /5 2 c 2 , Eq. (16) reduces to, 

r ~ — !— ■ . (17) 
1 2 — 1 

For r 2 = 4/3, r ~ 3, the standard ultrarelativistic result. 

To obtain a unique shock solution for 70 > 3, we modify both the shock structure and the 
compression ratio, rjjc, for the Monte Carlo simulation, calculate P 2 and e 2 from the resultant 
downstream particle distributions, and check that r^c — r as determined from Eq. (16). If Vmc r i 
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Fig. 12. — Particle distributions, p 4 f(p), for the shocks shown in Figs. 10 and 11 with 70 = 10. The 
spectra for the nonlinear (NL) and unmodified (UM) shocks are labeled and both are calculated at 
x = in the shock frame. The light-weight solid lines show spectral indexes, a = 4.23. 
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the shock structure and ryic are varied until a consistent solution is found with tmc — r - For smaller 
7o's, Eq. (16) doesn't apply because escaping fluxes become significant. In this case, however, the 
changes in the momentum and energy fluxes from changes in the shock structure and tmc are large 
enough that a unique solution can be found easily, as in our 70 = 1.4 example. 

Despite the fact that u(x) is modified on a fairly small length scale, the resultant particle 
distribution function is changed substantially, as indicated in Fig. 12. In this figure, the solid curve is 
from the unmodified shock and the dashed curve is from the nonlinear shock, both having exactly the 
same input conditions. The unmodified spectrum in Fig. 12 is similar to that shown with a dashed 
curve in Fig. 3 only now the far upstream plasma is taken to be a thermal gas at a temperature of 
10 6 K rather than a delta function distribution of particles with speeds, v m j = (2E- m j/m p ) 1 ^ 2 , with 
-Einj = 1 MeV, as was assumed for the example in Fig. 3. The shock smoothing has caused the low 
energy portion of the distribution to steepen and the overall intensity to decrease, reflecting the fact 
that the NL spectrum conserves momentum and energy while the UM one doesn't. The peaks in 
the two distributions at low momenta are also very different, with the NL spectrum having a larger 
fraction of slower particles than the UM one. These peaks result from the first shock crossing where 
all particles receive a large energy gain. In the UM case, a far greater fraction of the accelerated 
downstream particles are able to receive further energization by recrossing back into the upstream 
region than in the NL shock. The different speed distributions result in different T's and we find, 
by directly calculating T from the distributions, that Tum — 1-34 while Tnl — 1-36. 



4.3. Acceleration efficiency 

The absolute acceleration efficiency can be determined in our self-consistent, nonlinear examples 
directly from the particle distributions. In Fig. 13 we plot, e(> p), i.e., the fraction of kinetic energy 
density above a momentum p versus p for our 70 = 1.4 and 70 = 10 examples. The fraction of kinetic 
energy in the quasi-thermal part of the distribution can be determined from the relatively sharp 
fall off of the distributions at low momenta. This occurs near m p c for 70 = 1.4 and near 10 m p c for 
70 = 10. If we somewhat arbitrarily define the acceleration efficiency for these two examples to be 
€ia(> rripc) and eio(> 10m p c), respectively, we have €\a{> m p c) ~ 0.7 and eio(> Wm p c) ~ 0.6. 
Of course, the behavior of e(> p) depends on the particle spectrum from which it is derived and 
thus e(> p) is flatter for the 70 = 1.4 shock, with r ~ 4.5, than for the 70 = 10 shock, where 
r ~ 3.03 ±0.01. Furthermore, e(> p) depends strongly on p max i n the 70 = 1.4 shock where particle 
escape plays an important role in determining r. The maximum momentum has little influence 
for 70 = 10 because of the steep spectrum. The fraction of kinetic energy density above W 3 m p c is 
~ 0.3 for 70 = 1.4 and ~ 0.1 for 70 = 10. 
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Fig. 13. — Acceleration efficiency, e(> p), defined as the fraction of total kinetic energy density 
above p versus p. The sharp drop off at low momentum indicates the extent of the 'thermal' peak. 
The fraction of energy density above this drop off is a self-consistent measure of the production 
efficiency for 'superthermak particles. 
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5. Summary and Conclusions 

Particles gain energy in collisionless shocks by scattering nearly elastically off magnetic turbu- 
lence, back and forth between the converging plasmas upstream and downstream from the shock. 
While this basic shock acceleration physics is independent of the speed of the shock, the mathe- 
matical modeling of the process depends critically on whether or not the acceleration is efficient 
and whether or not particle speeds, v, are large compared to the shock speed, uq. The Monte Carlo 
techniques discussed here, which do not require v ^> uq, are well suited for the study of relativistic 
shocks, and for any shock where nonlinear effects are important and the energetic particles originate 
as thermal particles in the unshocked plasma. Except for computational limits, these techniques 
allow calculations of efficient particle acceleration in shocks of any Lorentz factor. 

As a check of our code, we have demonstrated that we obtain the well known, test-particle 
power laws in fully nonrelativistic and ultrarelativistic parallel shocks (Figs. 1 and 3). In trans- 
relativistic shocks, however, no such canonical results exist because the shock compression ratio, 
r, depends on the upstream conditions in a nonlinear fashion, even for test-particle shocks where 
all effects of accelerated particles are ignored (e.g., Ellison & Reynolds 1991). We show how r, 
and the resulting power law spectral index, a, vary through the trans-relativistic regime in Fig. 4, 
where r has been determined by balancing the momentum and energy fluxes across the shock in an 
iterative process. These test-particle results with self-consistently determined compression ratios 
and adiabatic indexes are in close agreement with the recent analytic results of Kirk et al. (2000) 
and Gallant (2002). 

Despite the fact that relativistic shock theory has concentrated almost exclusively on test- 
particle acceleration, it is likely that relativistic shocks are not test-particle but inject and accelerate 
particles efficiently. The reason is that regardless of the ambient far upstream conditions, all 
particles that are over taken by an ultrarelativistic shock will receive a large boosts in energy 
AE > 7omc 2 in their first shock crossing. Thus, virtually all of the particles in the downstream 
region of an unmodified shock are strongly relativistic with v ~ c. The ability to overtake the 
shock from downstream and be further accelerated depends only on the particle speed (Eq. 7) and 
the presence of magnetic waves or turbulence with sufficient power in wavelengths on the order of 
the particle gyroradii to isotropize the downstream distributions. It is generally assumed that the 
necessary magnetic turbulence is self-generated and if enough turbulence is generated to scatter 
high momentum particles (with very low densities) that constitute a test-particle power law, there 
should be enough generated to isotropize lower momentum particles (which carry the bulk of the 
density). If acceleration can occur at all, we believe it is likely to occur efficiently making it 
necessary to calculate the shock structure and particle acceleration self-consistently. Furthermore, 
if relativistic shock theory is to be applied to gamma-ray bursts, where high conversion efficiencies 
are generally assumed, nonlinear effects must be calculated. 

When energetic particles are generated in sufficient numbers, the conservation of momentum 
and energy requires that their backpressure modify the shock structure. Two basic effects occur: 



-27- 



a precursor is formed when the upstream plasma is slowed by the backpressure of the accelerated 
particles and the overall compression ratio changes from the test-particle value as a result of high 
energy particles escaping and/or a change in the shocked plasma's adiabatic index T. As indicated 
by our 70 = 1.4 example (Section 4.1), mildly relativistic shocks act as nonrelativistic ones showing 
a dramatic weakening of the subshock combined with a large increase in r (Fig. 8). These changes 
result in a particle distribution which is both steeper than the test-particle power law at low 
momenta and flatter at high momenta (Fig. 9). 

In faster shocks (i.e., 70 > 3), the initial test-particle spectrum is steep enough that particle 
escape is unimportant so only changes in T determine r (Eq. 16). In contrast to nonrelativistic 
shocks where the production of relativistic particles causes the compression ratio to increase, we 
show that r decreases smoothly to 3 as 70 increases and the fraction of fully relativistic shocked 
particles approaches one (Fig. 4). 

Our most important result is that efficient, mildly relativistic shocks do not produce particle 
spectra close to the so-called 'universal' power law having a ~ 4.3. This may be important for the 
interpretation of gamma-ray bursts since the internal shocks assumed responsible for converting 
the bulk kinetic energy of the fireball into internal particle energy may be mildly relativistic and 
the external shocks, believed responsible for producing gamma-ray burst afterglows, will always 
go through a mildly relativistic phase (see Piran 1999, for a comprehensive review of gamma-ray 
bursts). 

The authors thank M. G. Baring, L. O'C. Drury, F. C. Jones, and E. Parizot for helpful 
discussions. D.C.E. is grateful to the Department of Physics and Astronomy at Rice University 
and to the Dublin Institute for Advanced Studies for hosting visits where part of this work was 
done. 
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